1. Introduction¶
Maverik is currently facing a significant challenge in preparing comprehensive ROI documents for their new store openings. This challenge stems from the difficulty of accurately predicting first-year sales, which in turn hinders the company's strategic financial planning and decision-making processes, ultimately affecting their sustainable growth. To address this issue, the primary analytical goal is to develop a predictive model capable of providing precise daily forecasts for merchandise sales, food sales, unleaded fuel sales, and diesel fuel sales for new Maverik stores during their first year of operation.¶
The Modeling is dedicated to achieving this objective. Its purpose is to build a model that can forecast merchandise sales, food sales, unleaded fuel sales, and diesel fuel sales for the next 365 days using historical data as a basis. The notebook's focus is on evaluating various time series forecasting models to determine the most suitable option that aligns with Maverik's specific needs. Since there are four target variables involved, this constitutes a multivariate time series analysis. The two models tested are Prophet and XGBoost.¶
2. Exploratory Data Analysis¶
In this section, we provide an overview of the Exploratory Data Analysis carried out as an integral part of the project. Following that, we delve into the Modeling phase, where we assess different models to identify the most suitable one for our future predictions.¶
This section also encompasses strategies for data preparation, which involves managing variable transformations, conducting feature engineering, and addressing missing values in the dataset.¶
2.1. Importing Packages¶
In [74]:
#pip install xgboost
In [75]:
import numpy as np
import pandas as pd
import seaborn as sns
import warnings
warnings.filterwarnings("ignore")
import matplotlib.pyplot as plt
from pandas.plotting import table
import plotly.express as px
import matplotlib.ticker as mtick
In [76]:
import plotly.io as pio
pio.renderers.default = 'notebook_connected'
In [77]:
import xgboost as xgb
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import TimeSeriesSplit
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
In [78]:
from prophet import Prophet
2.2. Reading Data¶
In [79]:
q_data= pd.read_csv("qualitative_data_msba.csv")
t_data=pd.read_csv("time_series_data_msba.csv")
3. Summarizing & Understading the Data¶
In [80]:
q_data.head()
Out[80]:
In [81]:
t_data.head()
Out[81]:
3.1 Removing 1st column from both the dataframes named 'Unnamed: 0'¶
In [82]:
q_data= q_data.drop("Unnamed: 0", axis=1)
t_data= t_data.drop("Unnamed: 0", axis=1)
In [83]:
q_data.head()
Out[83]:
In [84]:
t_data.head()
Out[84]:
3.2 Checking the shape of both the dataframes¶
In [85]:
q_data.shape
Out[85]:
In [86]:
t_data.shape
Out[86]:
3.3 Obtaining description of the data in the datafame¶
In [87]:
q_data[['square_feet','x1_mile_pop','x1_mile_emp','x1_mile_income','x1_2_mile_pop','x1_2_mile_emp']].describe()
Out[87]:
In [88]:
t_data[['daily_yoy_ndt.total_inside_sales', 'daily_yoy_ndt.total_food_service', 'diesel', 'unleaded']].describe()
Out[88]:
3.4 Identified Zero Variance Variables¶
The following columns exhibit consistent values throughout and do not demonstrate significant correlation with the target variables. Therefore, they can be considered as zero variance variables and may be removed from the dataset:¶
1. front_door_count¶
2. godfather_s_pizza¶
3. diesel¶
4. car_wash¶
5. ev_charging¶
6. non_24_hour¶
7. self_check_out¶
In [89]:
q_data.drop(['front_door_count','godfather_s_pizza','car_wash','ev_charging','non_24_hour','self_check_out','diesel'], axis=1,inplace=True)
3.5 Checking the shape of q_data after columns drop¶
In [90]:
q_data.shape
Out[90]:
3.6 Finding NA and Categorical Values in Qualitative Dataframe (q_data)¶
In [91]:
na_values = q_data.isna().sum()
na_values[na_values > 0]
Out[91]:
As per the above result, 6 columns of qualitative dataframe has N/A value. It means that Maverik does not provide these features in few of it stores.¶
In [92]:
cat_col= q_data.select_dtypes(include="object").columns.tolist()
cat_col
Out[92]:
3.6.1 Replacing NA values with None¶
In [93]:
columns_with_na = q_data.columns[q_data.isna().any()].tolist()
for column in columns_with_na:
mode_value = q_data[column].mode()[0]
q_data[column].fillna("None", inplace=True)
3.6.2 Factorizing the Categorical columns¶
In [94]:
for column in cat_col:
q_data[column] = pd.factorize(q_data[column])[0]
In [95]:
import pandas as pd
# Define bins and labels
bins = [20, 30, 40, 50]
labels = [1, 2, 3]
# Bin the values in the 'parking_spaces' column
q_data['parking_spaces'] = pd.cut(q_data['parking_spaces'], bins=bins, labels=labels, include_lowest=True).astype(int)
In [96]:
q_data.head()
Out[96]:
3.7 Finding NA and Categorical Values in Time Series Dataframe (t_data)¶
In [97]:
na_values = t_data.isna().sum()
na_values[na_values > 0]
Out[97]:
There are no columns in time series data that has N/A value¶
3.8 Extracting Year, Month and Date from Calendar_day_date¶
In [98]:
t_data['calendar.calendar_day_date'] = pd.to_datetime(t_data['calendar.calendar_day_date'])
t_data['calendar_year'] = t_data['calendar.calendar_day_date'].dt.year
t_data['calendar_month'] = t_data['calendar.calendar_day_date'].dt.month
t_data['calendar_day'] = t_data['calendar.calendar_day_date'].dt.day
In [99]:
t_data.head()
Out[99]:
4. Analysis On Qualitative and Time Series Data¶
4.1 Analysing the 4 target variables(In-Store, Food, Diesel, Unleaded)¶
4.1.1 Plotting diesel and unleaded sales on Days of the week for the year 2021,2022 and 2023¶
In [100]:
t_data['calendar.calendar_day_date'] = pd.to_datetime(t_data['calendar.calendar_day_date'])
t_data.set_index('calendar.calendar_day_date', inplace=True)
Use the slider to drill through year,month,date and time¶
In [101]:
average_diesel_sales = t_data.groupby(t_data.index)['diesel'].mean()
average_unleaded_sales = t_data.groupby(t_data.index)['unleaded'].mean()
average_sales_df = pd.DataFrame({'Diesel': average_diesel_sales, 'Unleaded': average_unleaded_sales})
fig = px.line(average_sales_df, x=average_sales_df.index, y=average_sales_df.columns,
title='Average Diesel and Unleaded Sales Over Time')
fig.update_xaxes(title='Date', rangeslider_visible=True)
fig.update_yaxes(title='Average Sales')
fig.show()
4.1.2 Plotting Inside and Food sales on Days of the week for the year 2021,2022 and 2023¶
In [102]:
average_inside_sales = t_data.groupby(t_data.index)['daily_yoy_ndt.total_inside_sales'].mean()
average_food_sales = t_data.groupby(t_data.index)['daily_yoy_ndt.total_food_service'].mean()
average_sales_df = pd.DataFrame({'daily_yoy_ndt.total_inside_sales': average_inside_sales, 'daily_yoy_ndt.total_food_service': average_food_sales})
fig = px.line(average_sales_df, x=average_sales_df.index, y=average_sales_df.columns,
title='Average Inside and Food Sales Over Time')
fig.update_xaxes(title='Date', rangeslider_visible=True)
fig.update_yaxes(title='Average Sales')
fig.show()
4.1.3 Sales of all target variables based on Calender Week¶
In [103]:
t_data['calendar.fiscal_week_id_for_year'].min(), t_data['calendar.fiscal_week_id_for_year'].max()
Out[103]:
In [104]:
df_time_aggregate = t_data.groupby('calendar.fiscal_week_id_for_year')['daily_yoy_ndt.total_inside_sales', 'daily_yoy_ndt.total_food_service', 'diesel', 'unleaded'].sum().reset_index()
Lets check if we have more than one sales data for a date¶
In [105]:
temp = t_data.groupby('calendar.fiscal_week_id_for_year')['daily_yoy_ndt.total_inside_sales'].size()
temp[temp > 1].sort_values(ascending=False)
Out[105]:
Indexing with Time Series Data¶
In [106]:
df_time_aggregate = df_time_aggregate.set_index('calendar.fiscal_week_id_for_year')
df_time_aggregate.index
Out[106]:
In [107]:
plt.figure(figsize=(40, 50));
df_time_aggregate.plot();
4.1.4 Identifying Sale of Diesel and Unleaded during Weekend and Weekdays of the years 2021, 2022, 2023¶
In [108]:
for year in [2021, 2022, 2023]:
sales_year = t_data[t_data['calendar_year'] == year]
pivot_table = sales_year.pivot_table(index='calendar_information.type_of_day',
values=['diesel', 'unleaded'],
aggfunc='sum')
ax = pivot_table.plot(kind='bar', figsize=(10, 6))
plt.title(f'Diesel vs. Unleaded Sales by Type of Day - {year}')
plt.xlabel('Type of Day')
plt.ylabel('Total Sales')
ax.yaxis.set_major_formatter(mtick.FuncFormatter(lambda x, p: format(int(x), ',')))
ax.set_xticklabels(pivot_table.index, rotation=0)
plt.legend(['Diesel', 'Unleaded'])
plt.show()
4.1.5 Identifying Sale of Merchandise and Food during Weekend and Weekdays for the years 2021,2022,2023¶
In [109]:
for year in [2021, 2022, 2023]:
sales_year = t_data[t_data['calendar_year'] == year]
pivot_table = sales_year.pivot_table(index='calendar_information.type_of_day',
values=['daily_yoy_ndt.total_inside_sales', 'daily_yoy_ndt.total_food_service'],
aggfunc='sum')
ax = pivot_table.plot(kind='bar', figsize=(10, 6))
plt.title(f'Merchandise vs. Food Sales by Type of Day - {year}')
plt.xlabel('Type of Day')
plt.ylabel('Total Sales')
ax.yaxis.set_major_formatter(mtick.FuncFormatter(lambda x, p: format(int(x), ',')))
4.1.6 Plotting diesel and unleaded sales on Holidays for the year 2021,2022 and 2023¶
In [110]:
#2021
t_data_2021 = t_data[t_data['calendar_year'] == 2021]
sales_by_holiday_2021 = t_data_2021.groupby('calendar_information.holiday')[['diesel', 'unleaded']].sum()
sales_by_holiday_2021 = sales_by_holiday_2021[sales_by_holiday_2021.index != 'NONE']
holidays_2021 = sales_by_holiday_2021.index.tolist()
diesel_sales_2021 = sales_by_holiday_2021['diesel'].tolist()
unleaded_sales_2021 = sales_by_holiday_2021['unleaded'].tolist()
plt.figure(figsize=(15, 6))
bar_width = 0.35
r1 = range(len(holidays_2021))
r2 = [x + bar_width for x in r1]
plt.bar(r1, diesel_sales_2021, color='blue', width=bar_width, edgecolor='grey', label='Diesel')
plt.bar(r2, unleaded_sales_2021, color='orange', width=bar_width, edgecolor='grey', label='Unleaded')
plt.xlabel('Holiday')
plt.ylabel('Total Sales')
plt.title('Total Diesel and Unleaded Sales on Each Holiday in 2021')
plt.xticks([r + bar_width/2 for r in r1], holidays_2021, rotation=90)
plt.legend()
plt.show()
#2022
t_data_2022 = t_data[t_data['calendar_year'] == 2022]
sales_by_holiday_2022 = t_data_2022.groupby('calendar_information.holiday')[['diesel', 'unleaded']].sum()
sales_by_holiday_2022 = sales_by_holiday_2022[sales_by_holiday_2022.index != 'NONE']
holidays_2022 = sales_by_holiday_2022.index.tolist()
diesel_sales_2022 = sales_by_holiday_2022['diesel'].tolist()
unleaded_sales_2022 = sales_by_holiday_2022['unleaded'].tolist()
plt.figure(figsize=(15, 6))
bar_width = 0.35
r1 = range(len(holidays_2022))
r2 = [x + bar_width for x in r1]
plt.bar(r1, diesel_sales_2022, color='blue', width=bar_width, edgecolor='grey', label='Diesel')
plt.bar(r2, unleaded_sales_2022, color='orange', width=bar_width, edgecolor='grey', label='Unleaded')
plt.xlabel('Holiday')
plt.ylabel('Total Sales')
plt.title('Total Diesel and Unleaded Sales on Each Holiday in 2022')
plt.xticks([r + bar_width/2 for r in r1], holidays_2022, rotation=90)
plt.legend()
plt.show()
#2023
t_data_2023 = t_data[t_data['calendar_year'] == 2023]
sales_by_holiday_2023 = t_data_2023.groupby('calendar_information.holiday')[['diesel', 'unleaded']].sum()
sales_by_holiday_2023 = sales_by_holiday_2023[sales_by_holiday_2023.index != 'NONE']
holidays_2023 = sales_by_holiday_2023.index.tolist()
diesel_sales_2023 = sales_by_holiday_2023['diesel'].tolist()
unleaded_sales_2023 = sales_by_holiday_2023['unleaded'].tolist()
plt.figure(figsize=(15, 6))
bar_width = 0.35
r1 = range(len(holidays_2023))
r2 = [x + bar_width for x in r1]
plt.bar(r1, diesel_sales_2023, color='blue', width=bar_width, edgecolor='grey', label='Diesel')
plt.bar(r2, unleaded_sales_2023, color='orange', width=bar_width, edgecolor='grey', label='Unleaded')
plt.xlabel('Holiday')
plt.ylabel('Total Sales')
plt.title('Total Diesel and Unleaded Sales on Each Holiday in 2023')
plt.xticks([r + bar_width/2 for r in r1], holidays_2023, rotation=90)
plt.legend()
plt.show()
The above graphs shows that the sale of Unleaded is more than diesel on the Holidays.¶
4.1.7 Plotting in store and food sales on Holidays for the year 2021,2022 and 2023¶
In [111]:
#2021
t_data_2021 = t_data[t_data['calendar_year'] == 2021]
sales_by_holiday_2021 = t_data_2021.groupby('calendar_information.holiday')[['daily_yoy_ndt.total_inside_sales', 'daily_yoy_ndt.total_food_service']].sum()
sales_by_holiday_2021 = sales_by_holiday_2021[sales_by_holiday_2021.index != 'NONE']
holidays_2021 = sales_by_holiday_2021.index.tolist()
inStore_sales_2021 = sales_by_holiday_2021['daily_yoy_ndt.total_inside_sales'].tolist()
Food_sales_2021 = sales_by_holiday_2021['daily_yoy_ndt.total_food_service'].tolist()
plt.figure(figsize=(15, 6))
bar_width = 0.35
r1 = range(len(holidays_2021))
r2 = [x + bar_width for x in r1]
plt.bar(r1, inStore_sales_2021, color='blue', width=bar_width, edgecolor='grey', label='InStore')
plt.bar(r2, Food_sales_2021, color='orange', width=bar_width, edgecolor='grey', label='Food')
plt.xlabel('Holiday')
plt.ylabel('Total Sales')
plt.title('Total InStore and Food Sales on Each Holiday in 2021')
plt.xticks([r + bar_width/2 for r in r1], holidays_2021, rotation=90)
plt.legend()
plt.show()
#2022
t_data_2022 = t_data[t_data['calendar_year'] == 2022]
sales_by_holiday_2022 = t_data_2022.groupby('calendar_information.holiday')[['daily_yoy_ndt.total_inside_sales', 'daily_yoy_ndt.total_food_service']].sum()
sales_by_holiday_2022 = sales_by_holiday_2022[sales_by_holiday_2022.index != 'NONE']
holidays_2022 = sales_by_holiday_2022.index.tolist()
inStore_sales_2022 = sales_by_holiday_2022['daily_yoy_ndt.total_inside_sales'].tolist()
Food_sales_2022 = sales_by_holiday_2022['daily_yoy_ndt.total_food_service'].tolist()
plt.figure(figsize=(15, 6))
bar_width = 0.35
r1 = range(len(holidays_2022))
r2 = [x + bar_width for x in r1]
plt.bar(r1, inStore_sales_2022, color='blue', width=bar_width, edgecolor='grey', label='InStore')
plt.bar(r2, Food_sales_2022, color='orange', width=bar_width, edgecolor='grey', label='Food')
plt.xlabel('Holiday')
plt.ylabel('Total Sales')
plt.title('Total InStore and Food Sales on Each Holiday in 2022')
plt.xticks([r + bar_width/2 for r in r1], holidays_2022, rotation=90)
plt.legend()
plt.show()
#2023
t_data_2023 = t_data[t_data['calendar_year'] == 2023]
sales_by_holiday_2023 = t_data_2023.groupby('calendar_information.holiday')[['daily_yoy_ndt.total_inside_sales', 'daily_yoy_ndt.total_food_service']].sum()
sales_by_holiday_2023 = sales_by_holiday_2023[sales_by_holiday_2023.index != 'NONE']
holidays_2023 = sales_by_holiday_2023.index.tolist()
inStore_sales_2023 = sales_by_holiday_2023['daily_yoy_ndt.total_inside_sales'].tolist()
Food_sales_2023 = sales_by_holiday_2023['daily_yoy_ndt.total_food_service'].tolist()
plt.figure(figsize=(15, 6))
bar_width = 0.35
r1 = range(len(holidays_2023))
r2 = [x + bar_width for x in r1]
plt.bar(r1, inStore_sales_2023, color='blue', width=bar_width, edgecolor='grey', label='InStore')
plt.bar(r2, Food_sales_2023, color='orange', width=bar_width, edgecolor='grey', label='Food')
plt.xlabel('Holiday')
plt.ylabel('Total Sales')
plt.title('Total InStore and Food Sales on Each Holiday in 2021')
plt.xticks([r + bar_width/2 for r in r1], holidays_2023, rotation=90)
plt.legend()
plt.show()
The above graphs shows that the sale of In-store merchandise is more than the sale of food during Holidays.¶
4.1.8 Plotting In Store and food sales of Maverik store during each month of the year 2021,2022 and 2023¶
In [112]:
average_sales_2021 = t_data_2021.groupby('calendar_month')['daily_yoy_ndt.total_inside_sales'].mean()
average_sales_2022 = t_data_2022.groupby('calendar_month')['daily_yoy_ndt.total_inside_sales'].mean()
average_sales_2023 = t_data_2023.groupby('calendar_month')['daily_yoy_ndt.total_inside_sales'].mean()
months = range(1,13)
months_2023 = range(1, 9)
plt.figure(figsize=(10, 6))
plt.plot(months, average_sales_2021, label='2021')
plt.plot(months, average_sales_2022, label='2022')
plt.plot(months_2023, average_sales_2023, label='2023')
plt.xlabel('Month')
plt.ylabel('Average Sales')
plt.title('Average Inside Sales per Month (2021-2023)')
plt.legend()
plt.show()
average_sales_2021 = t_data_2021.groupby('calendar_month')['daily_yoy_ndt.total_food_service'].mean()
average_sales_2022 = t_data_2022.groupby('calendar_month')['daily_yoy_ndt.total_food_service'].mean()
average_sales_2023 = t_data_2023.groupby('calendar_month')['daily_yoy_ndt.total_food_service'].mean()
months_2023 = range(1, 9)
plt.figure(figsize=(10, 6))
plt.plot(months, average_sales_2021, label='2021')
plt.plot(months, average_sales_2022, label='2022')
plt.plot(months_2023, average_sales_2023, label='2023')
plt.xlabel('Month')
plt.ylabel('Average Sales')
plt.title('Average Food Sales per Month (2021-2023)')
plt.legend()
plt.show()
From the graphs depicted above, we can observe that the sale of in-store and food items increase during the end of Spring and whole summer season, it starts dipping during the end of Fall with the minimum sales observed during the last month of the year. We can roughly predict that this year sales are going to dip after Fall.¶
4.1.9 Plotting Diesel and Unleaded sales of Maverik store during each month of the year 2021,2022 and 2023¶
In [113]:
average_sales_2021 = t_data_2021.groupby('calendar_month')['diesel'].mean()
average_sales_2022 = t_data_2022.groupby('calendar_month')['diesel'].mean()
average_sales_2023 = t_data_2023.groupby('calendar_month')['diesel'].mean()
months_2023 = range(1, 9)
plt.figure(figsize=(10, 6))
plt.plot(months, average_sales_2021, label='2021')
plt.plot(months, average_sales_2022, label='2022')
plt.plot(months_2023, average_sales_2023, label='2023')
plt.xlabel('Month')
plt.ylabel('Average Sales')
plt.title('Average Diesel Sales per Month (2021-2023)')
plt.legend()
plt.show()
average_sales_2021 = t_data_2021.groupby('calendar_month')['unleaded'].mean()
average_sales_2022 = t_data_2022.groupby('calendar_month')['unleaded'].mean()
average_sales_2023 = t_data_2023.groupby('calendar_month')['unleaded'].mean()
months_2023 = range(1, 9)
plt.figure(figsize=(10, 6))
plt.plot(months, average_sales_2021, label='2021')
plt.plot(months, average_sales_2022, label='2022')
plt.plot(months_2023, average_sales_2023, label='2023')
plt.xlabel('Month')
plt.ylabel('Average Sales')
plt.title('Average unleaded Sales per Month (2021-2023)')
plt.legend()
plt.show()
We can see that year 2023 has the higher average sales of unleaded during July and August compared to the previous years whereas Diesel observed a good rise till July 2023 but there was a sudden dip in its average sales during August 2023.¶
4.2.0 Inside and Food Sales based on Seasonality¶
In [114]:
def get_season(month):
if month in [12, 1, 2]:
return 'Winter'
elif month in [3, 4, 5]:
return 'Spring'
elif month in [6, 7, 8]:
return 'Summer'
else:
return 'Fall'
t_data['season'] = t_data['calendar_month'].apply(get_season)
t_data_2021 = t_data[t_data['calendar_year'] == 2021].copy()
t_data_2022 = t_data[t_data['calendar_year'] == 2022].copy()
t_data_2023 = t_data[t_data['calendar_year'] == 2023].copy()
sales_by_season_2021 = t_data_2021.groupby('season')[['daily_yoy_ndt.total_inside_sales']].sum()
sales_by_season_2022 = t_data_2022.groupby('season')[['daily_yoy_ndt.total_inside_sales']].sum()
sales_by_season_2023 = t_data_2023.groupby('season')[['daily_yoy_ndt.total_inside_sales']].sum()
plt.figure(figsize=(18, 5))
plt.subplot(1, 3, 1)
plt.pie(sales_by_season_2021['daily_yoy_ndt.total_inside_sales'], labels=sales_by_season_2021.index, autopct='%1.1f%%', startangle=140, colors=['#ff9999','#66b3ff','#99ff99','#ffcc99'])
plt.title('Distribution of Inside Sales by Season in 2021')
plt.subplot(1, 3, 2)
plt.pie(sales_by_season_2022['daily_yoy_ndt.total_inside_sales'], labels=sales_by_season_2022.index, autopct='%1.1f%%', startangle=140, colors=['#ff9999','#66b3ff','#99ff99','#ffcc99'])
plt.title('Distribution of Inside Sales by Season in 2022')
plt.subplot(1, 3, 3)
plt.pie(sales_by_season_2023['daily_yoy_ndt.total_inside_sales'], labels=sales_by_season_2023.index, autopct='%1.1f%%', startangle=140, colors=['#66b3ff','#99ff99','#ffcc99'])
plt.title('Distribution of Inside Sales by Season in 2023')
plt.tight_layout()
sales_table_2021 = sales_by_season_2021.reset_index()
sales_table_2021.columns = ['Season', 'Total Sales (2021)']
sales_table_2022 = sales_by_season_2022.reset_index()
sales_table_2022.columns = ['Season', 'Total Sales (2022)']
sales_table_2023 = sales_by_season_2023.reset_index()
sales_table_2023.columns = ['Season', 'Total Sales (2023)']
with pd.option_context('display.float_format', '{:,.2f}'.format):
print("Inside Sales for 2021:")
display(sales_table_2021)
print("\n Inside Sales for 2022:")
display(sales_table_2022)
print("\n Inside Sales for 2023:")
display(sales_table_2023)
t_data['season'] = t_data['calendar_month'].apply(get_season)
t_data_2021 = t_data[t_data['calendar_year'] == 2021].copy()
t_data_2022 = t_data[t_data['calendar_year'] == 2022].copy()
t_data_2023 = t_data[t_data['calendar_year'] == 2023].copy()
sales_by_season_2021 = t_data_2021.groupby('season')[['daily_yoy_ndt.total_food_service']].sum()
sales_by_season_2022 = t_data_2022.groupby('season')[['daily_yoy_ndt.total_food_service']].sum()
sales_by_season_2023 = t_data_2023.groupby('season')[['daily_yoy_ndt.total_food_service']].sum()
plt.figure(figsize=(18, 5))
plt.subplot(1, 3, 1)
plt.pie(sales_by_season_2021['daily_yoy_ndt.total_food_service'], labels=sales_by_season_2021.index, autopct='%1.1f%%', startangle=140, colors=['#ff9999','#66b3ff','#99ff99','#ffcc99'])
plt.title('Distribution of Food Sales by Season in 2021')
plt.subplot(1, 3, 2)
plt.pie(sales_by_season_2022['daily_yoy_ndt.total_food_service'], labels=sales_by_season_2022.index, autopct='%1.1f%%', startangle=140, colors=['#ff9999','#66b3ff','#99ff99','#ffcc99'])
plt.title('Distribution of Food Sales by Season in 2022')
plt.subplot(1, 3, 3)
plt.pie(sales_by_season_2023['daily_yoy_ndt.total_food_service'], labels=sales_by_season_2023.index, autopct='%1.1f%%', startangle=140, colors=['#66b3ff','#99ff99','#ffcc99'])
plt.title('Distribution of Food Sales by Season in 2023')
plt.tight_layout()
sales_table_2021 = sales_by_season_2021.reset_index()
sales_table_2021.columns = ['Season', 'Total Sales (2021)']
sales_table_2022 = sales_by_season_2022.reset_index()
sales_table_2022.columns = ['Season', 'Total Sales (2022)']
sales_table_2023 = sales_by_season_2023.reset_index()
sales_table_2023.columns = ['Season', 'Total Sales (2023)']
with pd.option_context('display.float_format', '{:,.2f}'.format):
print("Food Sales for 2021:")
display(sales_table_2021)
print("\n Food Sales for 2022:")
display(sales_table_2022)
print("\n Food Sales for 2023:")
display(sales_table_2023)
4.2.1 Diesel and Unleaded Sales based on Seasonality¶
In [115]:
def get_season(month):
if month in [12, 1, 2]:
return 'Winter'
elif month in [3, 4, 5]:
return 'Spring'
elif month in [6, 7, 8]:
return 'Summer'
else:
return 'Fall'
t_data['season'] = t_data['calendar_month'].apply(get_season)
t_data_2021 = t_data[t_data['calendar_year'] == 2021].copy()
t_data_2022 = t_data[t_data['calendar_year'] == 2022].copy()
t_data_2023 = t_data[t_data['calendar_year'] == 2023].copy()
sales_by_season_2021 = t_data_2021.groupby('season')[['diesel']].sum()
sales_by_season_2022 = t_data_2022.groupby('season')[['diesel']].sum()
sales_by_season_2023 = t_data_2023.groupby('season')[['diesel']].sum()
plt.figure(figsize=(18, 5))
plt.subplot(1, 3, 1)
plt.pie(sales_by_season_2021['diesel'], labels=sales_by_season_2021.index, autopct='%1.1f%%', startangle=140, colors=['#ff9999','#66b3ff','#99ff99','#ffcc99'])
plt.title('Distribution of Sales by Season in 2021')
plt.subplot(1, 3, 2)
plt.pie(sales_by_season_2022['diesel'], labels=sales_by_season_2022.index, autopct='%1.1f%%', startangle=140, colors=['#ff9999','#66b3ff','#99ff99','#ffcc99'])
plt.title('Distribution of Sales by Season in 2022')
plt.subplot(1, 3, 3)
plt.pie(sales_by_season_2023['diesel'], labels=sales_by_season_2023.index, autopct='%1.1f%%', startangle=140, colors=['#66b3ff','#99ff99','#ffcc99'])
plt.title('Distribution of Sales by Season in 2023')
plt.tight_layout()
sales_table_2021 = sales_by_season_2021.reset_index()
sales_table_2021.columns = ['Season', 'Total Sales (2021)']
sales_table_2022 = sales_by_season_2022.reset_index()
sales_table_2022.columns = ['Season', 'Total Sales (2022)']
sales_table_2023 = sales_by_season_2023.reset_index()
sales_table_2023.columns = ['Season', 'Total Sales (2023)']
with pd.option_context('display.float_format', '{:,.2f}'.format):
print("Sales for 2021:")
display(sales_table_2021)
print("\nSales for 2022:")
display(sales_table_2022)
print("\nSales for 2023:")
display(sales_table_2023)
t_data['season'] = t_data['calendar_month'].apply(get_season)
t_data_2021 = t_data[t_data['calendar_year'] == 2021].copy()
t_data_2022 = t_data[t_data['calendar_year'] == 2022].copy()
t_data_2023 = t_data[t_data['calendar_year'] == 2023].copy()
sales_by_season_2021 = t_data_2021.groupby('season')[['unleaded']].sum()
sales_by_season_2022 = t_data_2022.groupby('season')[['unleaded']].sum()
sales_by_season_2023 = t_data_2023.groupby('season')[['unleaded']].sum()
plt.figure(figsize=(18, 5))
plt.subplot(1, 3, 1)
plt.pie(sales_by_season_2021['unleaded'], labels=sales_by_season_2021.index, autopct='%1.1f%%', startangle=140, colors=['#ff9999','#66b3ff','#99ff99','#ffcc99'])
plt.title('Distribution of Sales by Season in 2021')
plt.subplot(1, 3, 2)
plt.pie(sales_by_season_2022['unleaded'], labels=sales_by_season_2022.index, autopct='%1.1f%%', startangle=140, colors=['#ff9999','#66b3ff','#99ff99','#ffcc99'])
plt.title('Distribution of Sales by Season in 2022')
plt.subplot(1, 3, 3)
plt.pie(sales_by_season_2023['unleaded'], labels=sales_by_season_2023.index, autopct='%1.1f%%', startangle=140, colors=['#66b3ff','#99ff99','#ffcc99'])
plt.title('Distribution of Sales by Season in 2023')
plt.tight_layout()
sales_table_2021 = sales_by_season_2021.reset_index()
sales_table_2021.columns = ['Season', 'Total Sales (2021)']
sales_table_2022 = sales_by_season_2022.reset_index()
sales_table_2022.columns = ['Season', 'Total Sales (2022)']
sales_table_2023 = sales_by_season_2023.reset_index()
sales_table_2023.columns = ['Season', 'Total Sales (2023)']
with pd.option_context('display.float_format', '{:,.2f}'.format):
print("Sales for 2021:")
display(sales_table_2021)
print("\nSales for 2022:")
display(sales_table_2022)
print("\nSales for 2023:")
display(sales_table_2023)
4.2.2 Total sale of target variables for all site id's¶
In [116]:
grouped_data = t_data.groupby('site_id_msba').agg({
'daily_yoy_ndt.total_inside_sales': 'sum',
'daily_yoy_ndt.total_food_service': 'sum',
'diesel': 'sum',
'unleaded': 'sum'
}).reset_index()
grouped_data.set_index('site_id_msba', inplace=True)
ax = grouped_data.plot(kind='bar', figsize=(12, 8), width=0.8)
plt.xlabel('site_id_msba')
plt.ylabel('Total Sum')
plt.title('Total Sum of Sales and Fuel Types by site_id_msba')
plt.legend(title='Categories')
plt.tight_layout()
plt.show()
The above graph and table shows that site id 22085 recorded maximum sale of Merchandise and food service, site id 21980 recorded the maximum sale of Diesel and Site id 22260 recorded the maximum sale od Unleaded.¶
5. Merging both dataframes¶
In [117]:
df_merged=pd.merge(t_data,q_data, on='site_id_msba', how='left')
In [118]:
df_merged.shape
Out[118]:
In [119]:
df_merged.head()
Out[119]:
5.1 Relation between customer population and income¶
In [120]:
plt.scatter(df_merged['x1_mile_pop'], df_merged['x1_mile_income'])
plt.xlabel('Population within 1-mile radius')
plt.ylabel('Median income of 1-mile radius population')
plt.title('Scatter Plot: 1-mile Population vs. 1-mile Income')
plt.show()
plt.scatter(df_merged['x1_2_mile_pop'], df_merged['x1_2_mile_income'])
plt.xlabel('Population within 1/2-mile radius')
plt.ylabel('Median income of 1/2-mile radius population')
plt.title('Scatter Plot: 1/2-mile Population vs. 1/2-mile Income')
plt.show()
plt.scatter(df_merged['x5_min_pop'], df_merged['x5_min_inc'])
plt.xlabel('Population within 5-minute radius')
plt.ylabel('Median income of 5-minute radius population')
plt.title('Scatter Plot: 5-minute Population vs. 5-minute Income')
plt.show()
plt.scatter(df_merged['x7_min_pop'], df_merged['x7_min_inc'])
plt.xlabel('Population within 7-minute radius')
plt.ylabel('Median income of 7-minute radius population')
plt.title('Scatter Plot: 7-minute Population vs. 7-minute Income')
plt.show()
5.2 Population in Miles vs Total Sales¶
In [121]:
df_merged['total_sales'] = df_merged['daily_yoy_ndt.total_food_service'] + df_merged['daily_yoy_ndt.total_inside_sales']
total_sales_1_mile = df_merged[df_merged['x1_mile_pop'] > 0]['total_sales'].sum()
total_sales_2_mile = df_merged[df_merged['x1_2_mile_pop'] > 0]['total_sales'].sum()
total_sales_5_min = df_merged[df_merged['x5_min_pop'] > 0]['total_sales'].sum()
total_sales_7_min = df_merged[df_merged['x7_min_pop'] > 0]['total_sales'].sum()
labels = ['1 Mile', '2 Miles', '5 Minutes', '7 Minutes']
values = [total_sales_1_mile, total_sales_2_mile, total_sales_5_min, total_sales_7_min]
plt.figure(figsize=(8, 8))
plt.pie(values, labels=labels, autopct='%1.1f%%', startangle=140)
plt.title('Contribution of Total Sales by Population Radius')
plt.show()
5.3 Customer Income vs Total Sales¶
In [122]:
income_bins = [20000, 30000, 50000, 70000, 90000, 110000, 1000000]
income_labels = ['20k-30k', '30k-50k', '50k-70k', '70k-90k', '90k-110k', '100k+']
df_merged['x1_mile_income_bin'] = pd.cut(df_merged['x1_mile_income'], bins=income_bins, labels=income_labels, right=False)
df_merged['x1_2_mile_income_bin'] = pd.cut(df_merged['x1_2_mile_income'], bins=income_bins, labels=income_labels, right=False)
df_merged['x5_min_inc_bin'] = pd.cut(df_merged['x5_min_inc'], bins=income_bins, labels=income_labels, right=False)
df_merged['x7_min_inc_bin'] = pd.cut(df_merged['x7_min_inc'], bins=income_bins, labels=income_labels, right=False)
df_merged['total_sales'] = df_merged['daily_yoy_ndt.total_inside_sales'] + df_merged['daily_yoy_ndt.total_food_service']
total_sales_by_income_x1_mile = df_merged.groupby('x1_mile_income_bin')['total_sales'].sum()
total_sales_by_income_x1_2_mile = df_merged.groupby('x1_2_mile_income_bin')['total_sales'].sum()
total_sales_by_income_x5_min = df_merged.groupby('x5_min_inc_bin')['total_sales'].sum()
total_sales_by_income_x7_min = df_merged.groupby('x7_min_inc_bin')['total_sales'].sum()
plt.figure(figsize=(12, 8))
plt.subplot(221)
total_sales_by_income_x1_mile.plot(kind='bar')
plt.title('Total Sales by x1_mile_income_bin')
plt.xlabel('Income Range')
plt.ylabel('Total Sales')
plt.subplot(222)
total_sales_by_income_x1_2_mile.plot(kind='bar')
plt.title('Total Sales by x1_2_mile_income_bin')
plt.xlabel('Income Range')
plt.ylabel('Total Sales')
plt.subplot(223)
total_sales_by_income_x5_min.plot(kind='bar')
plt.title('Total Sales by x5_min_inc_bin')
plt.xlabel('Income Range')
plt.ylabel('Total Sales')
plt.subplot(224)
total_sales_by_income_x7_min.plot(kind='bar')
plt.title('Total Sales by x7_min_inc_bin')
plt.xlabel('Income Range')
plt.ylabel('Total Sales')
plt.tight_layout()
plt.show()
Customers with income range between 30k to 70k contribute majorly towards total sales¶
5.4 Distribution of Sales based on Stores that Sell "freal," "bonfire_grill," "pizza," and "cinnabon"¶
In [123]:
sells_items = df_merged[(df_merged['freal'] > 0) | (df_merged['bonfire_grill'] > 0) | (df_merged['pizza'] > 0) | (df_merged['cinnabon'] > 0)]
does_not_sell_items = df_merged[(df_merged['freal'] == 0) & (df_merged['bonfire_grill'] == 0) & (df_merged['pizza'] == 0) & (df_merged['cinnabon'] == 0)]
avg_food_service_sells_items = sells_items['daily_yoy_ndt.total_food_service'].mean()
avg_food_service_does_not_sell_items = does_not_sell_items['daily_yoy_ndt.total_food_service'].mean()
store_type_counts = [len(sells_items), len(does_not_sell_items)]
labels = ['Sells Items', 'Does Not Sell Items']
plt.figure(figsize=(8, 6))
plt.pie(store_type_counts, labels=labels, autopct='%1.1f%%', startangle=140)
plt.title('Distribution of Store Types')
plt.show()
We can see that 70.3% of the food sales come from stores that sell food items and 29.7% from those that do not sell food items¶
5.5 Comparing Sale of Diesel and Unleaded based on Fuel Station Features¶
In [124]:
columns_to_visualize = ['cat_scales', 'propane', 'traditional_forecourt_fueling_positions',
'traditional_forecourt_layout', 'traditional_forecourt_stack_type',
'hi_flow_rv_lanes_layout', 'hi_flow_rv_lanes_stack_type']
for column in columns_to_visualize:
sales_by_column = df_merged.groupby(column)[['diesel', 'unleaded']].sum()
sales_by_column.plot(kind='bar', figsize=(10, 6))
plt.title(f'Sales of Diesel and Unleaded by {column}')
plt.xlabel(column)
plt.ylabel('Total Sales')
plt.xticks(rotation=0)
plt.legend(['Diesel', 'Unleaded'])
plt.show()
In [125]:
### comment to be added by Shivi
5.6 Comparing Sale of Diesel based on Fuel Station Features¶
In [126]:
numerical_columns = ['hi_flow_lanes', 'hi_flow_rv_lanes', 'def',
'hi_flow_lanes_fueling_positions', 'hi_flow_lanes_layout',
'hi_flow_lanes_stack_type', 'diesel']
subset_df = df_merged[numerical_columns]
correlation_matrix = subset_df.corr()
plt.figure(figsize=(10, 8))
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', linewidths=0.5)
plt.title('Correlation Heatmap')
plt.show()
We observe that sale of diesel increases if there is an increase in the number of hi flow lanes fueling positions.¶
5.7 Comparing Sale of Unleaded based on Fuel Station Features¶
In [127]:
numerical_columns = [ 'ethanol_free', 'rv_lanes', 'rv_dumps',
'rv_lanes_fueling_positions', 'rv_lanes_layout',
'rv_lanes_stack_type', 'unleaded']
subset_df = df_merged[numerical_columns]
correlation_matrix = subset_df.corr()
plt.figure(figsize=(10, 8))
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', linewidths=0.5)
plt.title('Correlation Heatmap')
plt.show()
As reflected in the heat map, all the predictors except rv_dumps has a mild positive correlation where an increase in one unit of the value will lead to the increase in sale of unleaded¶
Dropping the columns not required during Modeling Phase¶
In [128]:
columns_to_drop = ['x1_mile_income_bin', 'x1_2_mile_income_bin', 'x5_min_inc_bin', 'x7_min_inc_bin']
df_merged = df_merged.drop(columns=columns_to_drop)
6.1 Feature Engineering¶
In [129]:
df_merged['calendar.calendar_day_date'] = t_data.index
In [130]:
df_merged.set_index('calendar.calendar_day_date', inplace=True)
One Hot Encoding¶
In [131]:
# Identifying categorical columns
categorical_columns = df_merged.select_dtypes(include=['object']).columns
# Applying one-hot encoding to categorical columns
encoded_df = pd.get_dummies(df_merged, columns=categorical_columns)
In [132]:
encoded_df.head()
Out[132]:
'encoded_df' now contains the original columns along with one-hot encoded categorical columns¶
3. Modeling process¶
I. XGBoost Model¶
A.Target Variable = Food¶
Plotting the graph for sale of food over the years¶
In [133]:
average_sales = encoded_df['daily_yoy_ndt.total_food_service'].mean()
encoded_df['daily_yoy_ndt.total_food_service'].plot(style=".", figsize=(15,5), title="Food Sale over the Years")
plt.axhline(y=average_sales, color='r', linestyle='-', label='Average Sales')
plt.legend()
plt.show()
Train Test Split¶
In [134]:
#splitting train as date less than 1st jan 2023 and test as later
train= encoded_df['daily_yoy_ndt.total_food_service'].loc[encoded_df.index< '2023-01-01']
test= encoded_df['daily_yoy_ndt.total_food_service'].loc[encoded_df.index>= '2023-01-01']
fig, ax= plt.subplots(figsize=(15,5))
train.plot(ax=ax, label='training set')
test.plot(ax=ax, label='Test set')
ax.axvline('2023-01-01',color='black',ls='--')
ax.legend(['Training Set','Test Set'])
plt.show()
Splitting data into training and testing sets¶
In [135]:
# Features (X) and Target (y)
X = encoded_df.drop(columns=['daily_yoy_ndt.total_food_service'])
y = encoded_df['daily_yoy_ndt.total_food_service']
train_end_date = '2023-01-01'
X_train = X[X.index < train_end_date]
X_test = X[X.index >= train_end_date]
y_train = y[y.index < train_end_date]
y_test = y[y.index >= train_end_date]
We now have X_train, X_test, y_train, and y_test ready for modeling¶
In [136]:
reg = xgb.XGBRegressor(base_score=0.5, booster='gbtree',
n_estimators=1000,
early_stopping_rounds=50,
objective='reg:linear',
max_depth=3,
learning_rate=0.01)
reg.fit(X_train, y_train,
eval_set=[(X_train, y_train), (X_test, y_test)],
verbose=100)
Out[136]:
Identifying Top 10 features impacting the Food Sale¶
In [137]:
fi = pd.DataFrame(data=reg.feature_importances_,
index=reg.feature_names_in_,
columns=['importance'])
top_10_features = fi.sort_values(by='importance', ascending=False).head(10)
top_10_features.plot(kind='barh', title='Top 10 Feature Importance')
Out[137]:
Plotting the graph to depict Train and Test Split¶
In [138]:
train = encoded_df[['daily_yoy_ndt.total_food_service','x1_2_mile_pop','x7_min_pop','x1_mile_pop','womens_sink_count','hi_flow_lanes','square_feet','x1_2_mile_emp','x5_min_pop']].loc[encoded_df.index < '2023-01-01']
test = encoded_df[['daily_yoy_ndt.total_food_service','x1_2_mile_pop','x7_min_pop','x1_mile_pop','womens_sink_count','hi_flow_lanes','square_feet','x1_2_mile_emp','x5_min_pop']].loc[encoded_df.index >= '2023-01-01']
fig, ax = plt.subplots(figsize=(15, 5))
train['daily_yoy_ndt.total_food_service'].plot(ax=ax, label='Training Set', title='Data Train/Test Split')
test['daily_yoy_ndt.total_food_service'].plot(ax=ax, label='Test Set')
ax.axvline('2023-01-01', color='black', ls='--')
ax.legend(['Training Set', 'Test Set'])
plt.show()
Time series cross validation¶
In [139]:
tss = TimeSeriesSplit(n_splits=5)
df = encoded_df[['daily_yoy_ndt.total_food_service','x1_2_mile_pop','x7_min_pop','x1_mile_pop','womens_sink_count','hi_flow_lanes','square_feet','x1_2_mile_emp','x5_min_pop']].sort_index()
In [140]:
fig, axs = plt.subplots(5, 1, figsize=(15, 15), sharex=True)
fold = 0
for train_idx, val_idx in tss.split(df):
train = df.iloc[train_idx]
test = df.iloc[val_idx]
train['daily_yoy_ndt.total_food_service'].plot(ax=axs[fold],
label='Training Set',
title=f'Data Train/Test Split Fold {fold}')
test['daily_yoy_ndt.total_food_service'].plot(ax=axs[fold],
label='Test Set')
axs[fold].axvline(test.index.min(), color='black', ls='--')
fold += 1
plt.show()
Creating time series features based on time series index¶
In [141]:
def create_features(df):
df = df.copy()
df['dayofweek'] = df.index.dayofweek
df['quarter'] = df.index.quarter
df['month'] = df.index.month
df['year'] = df.index.year
df['dayofyear'] = df.index.dayofyear
df['dayofmonth'] = df.index.day
return df
df = create_features(df)
In [142]:
df.head()
Out[142]:
Defining Lags¶
In [143]:
def add_lags(df):
target_map= df['daily_yoy_ndt.total_food_service'].to_dict()
df['lag1'] = (df.index - pd.Timedelta('364 days')).map(target_map)
return df
In [144]:
df = add_lags(df)
In [145]:
df.tail()
Out[145]:
Fitting the Model to calculate Evaluation Metrics¶
In [146]:
tss = TimeSeriesSplit(n_splits=5)
df = df.sort_index()
fold = 0
preds = []
scores = []
for train_idx, val_idx in tss.split(df):
train = df.iloc[train_idx]
test = df.iloc[val_idx]
TARGET = 'daily_yoy_ndt.total_food_service'
X_train = train
y_train = train[TARGET]
X_test = test
y_test = test[TARGET]
reg = xgb.XGBRegressor(base_score=0.5, booster='gbtree',
n_estimators=1000,
early_stopping_rounds=50,
objective='reg:linear',
max_depth=3,
learning_rate=0.01)
reg.fit(X_train, y_train,
eval_set=[(X_train, y_train), (X_test, y_test)],
verbose=100)
y_pred = reg.predict(X_test)
preds.append(y_pred)
score = np.sqrt(mean_squared_error(y_test, y_pred))
scores.append(score)
Calculation of evaluation metrics¶
In [147]:
mae = mean_absolute_error(y_test, y_pred)
mse = mean_squared_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)
# Print the evaluation metrics
print(f'Mean Absolute Error (MAE): {mae:.2f}')
print(f'Mean Squared Error (MSE): {mse:.2f}')
print(f'R-squared (R2) Score: {r2:.2f}')
In [148]:
print(f'Score across folds {np.mean(scores):0.4f}')
print(f'Fold scores:{scores}')
Retraining all data¶
In [149]:
future = pd.date_range('2023-08-17','2024-12-31', freq='D')
future_df = pd.DataFrame(index=future)
df['isFuture'] = False
TARGET = 'daily_yoy_ndt.total_food_service'
# Define features (excluding the target and 'isFuture' column)
FEATURES = ['x1_2_mile_pop','x7_min_pop','x1_mile_pop','womens_sink_count','hi_flow_lanes','square_feet','x1_2_mile_emp','x5_min_pop',
'dayofweek', 'quarter', 'month', 'year', 'dayofyear',
'dayofmonth', 'lag1','isFuture']
# Split the data into features (X) and target (y)
X = df[FEATURES]
y = df[TARGET]
# Create the XGBoost model
reg = xgb.XGBRegressor(base_score=0.5,
booster='gbtree',
n_estimators=500,
objective='reg:linear',
max_depth=3,
learning_rate=0.01)
# Fit the model
reg.fit(X, y, eval_set=[(X, y)], verbose=100)
Out[149]:
Creating Future Dataframe¶
In [150]:
future = pd.date_range('2023-08-17','2024-12-31', freq='D')
future_df = pd.DataFrame(index=future)
future_df['isFuture'] = True
df['isFuture'] = False
df_and_future = pd.concat([df, future_df])
df_and_future = create_features(df_and_future)
df_and_future = add_lags(df_and_future)
In [151]:
future_w_features = df_and_future.query('isFuture').copy()
In [152]:
future_w_features.head()
Out[152]:
In [153]:
future_w_features = future_w_features.drop('daily_yoy_ndt.total_food_service', axis=1)
In [154]:
future_w_features['pred'] = reg.predict(future_w_features)
In [155]:
future_w_features.index = pd.to_datetime(future_w_features.index) # Make sure index is datetime
plt.figure(figsize=(10, 5))
plt.plot(future_w_features.index, future_w_features['pred'], color='r', marker='o', linestyle='-', linewidth=1, markersize=3)
plt.title('Future Predictions For Food Sales ')
plt.xlabel('Date')
plt.ylabel('Inside Sales')
plt.grid(True)
plt.show()
Combining Future Prediction with the Historical Data¶
In [156]:
# Plotting historical data
plt.figure(figsize=(15, 5))
# Plot the entire historical data
encoded_df['daily_yoy_ndt.total_food_service'].plot(label='Historical Data', color='b')
plt.axvline('2023-08-17', color='black', ls='--')
# Plot future predictions
plt.plot(future_w_features.index, future_w_features['pred'], label='Future Predictions', color='red', linestyle='dashed')
# Set title and labels
plt.title('Historical Data and Future Predictions of Total Food Sales')
plt.xlabel('Date')
plt.ylabel('Value')
# Add legend
plt.legend()
# Show plot
plt.show()
Displaying the data for future predictions¶
In [157]:
start_date = '2023-08-17'
end_date = pd.to_datetime(start_date) + pd.DateOffset(days=502)
date_range = pd.date_range(start=start_date, end=end_date, freq='D')
# Create a DataFrame for the forecasted sales
forecasted_sales_df = pd.DataFrame({'Date': date_range, 'Sales': future_w_features['pred']})
# Display the DataFrame
print(forecasted_sales_df)
B.Target Variable = Inside Sales¶
Plotting the graph for sale of food over the years¶
In [158]:
average_sales = encoded_df['daily_yoy_ndt.total_inside_sales'].mean()
encoded_df['daily_yoy_ndt.total_inside_sales'].plot(style=".", figsize=(15,5), title="Diesel Sale over the Years")
plt.axhline(y=average_sales, color='r', linestyle='-', label='Average Sales')
plt.legend()
plt.show()
Train Test Split¶
In [159]:
train= encoded_df['daily_yoy_ndt.total_inside_sales'].loc[encoded_df.index< '2023-01-01']
test= encoded_df['daily_yoy_ndt.total_inside_sales'].loc[encoded_df.index>= '2023-01-01']
fig, ax= plt.subplots(figsize=(15,5))
train.plot(ax=ax, label='training set')
test.plot(ax=ax, label='Test set')
ax.axvline('2023-01-01',color='black',ls='--')
ax.legend(['Training Set','Test Set'])
plt.show()
Splitting data into training and testing sets¶
In [160]:
# Features (X) and Target (y)
X = encoded_df.drop(columns=['daily_yoy_ndt.total_inside_sales'])
y = encoded_df['daily_yoy_ndt.total_inside_sales']
# Split data into training and testing sets
train_end_date = '2023-01-01'
X_train = X[X.index < train_end_date]
X_test = X[X.index >= train_end_date]
y_train = y[y.index < train_end_date]
y_test = y[y.index >= train_end_date]
We now have X_train, X_test, y_train, and y_test ready for modeling¶
In [161]:
reg = xgb.XGBRegressor(base_score=0.5, booster='gbtree',
n_estimators=1000,
early_stopping_rounds=50,
objective='reg:linear',
max_depth=3,
learning_rate=0.01)
reg.fit(X_train, y_train,
eval_set=[(X_train, y_train), (X_test, y_test)],
verbose=100)
Out[161]:
Identifying Top 10 features impacting the Inside Sale¶
In [162]:
fi = pd.DataFrame(data=reg.feature_importances_,
index=reg.feature_names_in_,
columns=['importance'])
top_10_features = fi.sort_values(by='importance', ascending=False).head(10)
top_10_features.plot(kind='barh', title='Top 10 Feature Importance')
Out[162]:
Plotting the graph to depict Train and Test Split¶
In [163]:
train = encoded_df[['daily_yoy_ndt.total_inside_sales','x5_min_pop','x1_2_mile_emp','square_feet','hi_flow_lanes','womens_sink_count','x1_mile_pop','x7_min_pop','x1_2_mile_pop']].loc[encoded_df.index < '2023-01-01']
test = encoded_df[['daily_yoy_ndt.total_inside_sales','x5_min_pop','x1_2_mile_emp','square_feet','hi_flow_lanes','womens_sink_count','x1_mile_pop','x7_min_pop','x1_2_mile_pop']].loc[encoded_df.index >= '2023-01-01']
fig, ax = plt.subplots(figsize=(15, 5))
train['daily_yoy_ndt.total_inside_sales'].plot(ax=ax, label='Training Set', title='Data Train/Test Split')
test['daily_yoy_ndt.total_inside_sales'].plot(ax=ax, label='Test Set')
ax.axvline('2023-01-01', color='black', ls='--')
ax.legend(['Training Set', 'Test Set'])
plt.show()
Time series cross validation¶
In [164]:
tss = TimeSeriesSplit(n_splits=5)
df = encoded_df[['daily_yoy_ndt.total_inside_sales','x5_min_pop','x1_2_mile_emp','square_feet','hi_flow_lanes','womens_sink_count','x1_mile_pop','x7_min_pop','x1_2_mile_pop']].sort_index()
In [165]:
fig, axs = plt.subplots(5, 1, figsize=(15, 15), sharex=True)
fold = 0
for train_idx, val_idx in tss.split(df):
train = df.iloc[train_idx]
test = df.iloc[val_idx]
train['daily_yoy_ndt.total_inside_sales'].plot(ax=axs[fold],
label='Training Set',
title=f'Data Train/Test Split Fold {fold}')
test['daily_yoy_ndt.total_inside_sales'].plot(ax=axs[fold],
label='Test Set')
axs[fold].axvline(test.index.min(), color='black', ls='--')
fold += 1
plt.show()
Creating time series features based on time series index¶
In [166]:
def create_features(df):
df = df.copy()
df['dayofweek'] = df.index.dayofweek
df['quarter'] = df.index.quarter
df['month'] = df.index.month
df['year'] = df.index.year
df['dayofyear'] = df.index.dayofyear
df['dayofmonth'] = df.index.day
return df
df = create_features(df)
Defining Lags¶
In [167]:
def add_lags(df):
target_map= df['daily_yoy_ndt.total_inside_sales'].to_dict()
df['lag1'] = (df.index - pd.Timedelta('364 days')).map(target_map)
return df
In [168]:
df = add_lags(df)
In [169]:
df.tail()
Out[169]:
Fitting the Model to calculate the Evaluation Metrics¶
In [170]:
tss = TimeSeriesSplit(n_splits=5)
df = df.sort_index()
fold = 0
preds = []
scores = []
for train_idx, val_idx in tss.split(df):
train = df.iloc[train_idx]
test = df.iloc[val_idx]
TARGET = 'daily_yoy_ndt.total_inside_sales'
X_train = train
y_train = train[TARGET]
X_test = test
y_test = test[TARGET]
reg = xgb.XGBRegressor(base_score=0.5, booster='gbtree',
n_estimators=1000,
early_stopping_rounds=50,
objective='reg:linear',
max_depth=3,
learning_rate=0.01)
reg.fit(X_train, y_train,
eval_set=[(X_train, y_train), (X_test, y_test)],
verbose=100)
y_pred = reg.predict(X_test)
preds.append(y_pred)
score = np.sqrt(mean_squared_error(y_test, y_pred))
scores.append(score)
Calculation of evaluation metrics¶
In [171]:
mae = mean_absolute_error(y_test, y_pred)
mse = mean_squared_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)
# Print the evaluation metrics
print(f'Mean Absolute Error (MAE): {mae:.2f}')
print(f'Mean Squared Error (MSE): {mse:.2f}')
print(f'R-squared (R2) Score: {r2:.2f}')
In [172]:
print(f'Score across folds {np.mean(scores):0.4f}')
print(f'Fold scores:{scores}')
Retraining all data¶
In [173]:
future = pd.date_range('2023-08-17','2024-12-31', freq='D')
future_df = pd.DataFrame(index=future)
df['isFuture'] = False
TARGET = 'daily_yoy_ndt.total_inside_sales'
# Define features (excluding the target and 'isFuture' column)
FEATURES = ['x5_min_pop', 'x1_2_mile_emp', 'square_feet', 'hi_flow_lanes',
'womens_sink_count', 'x1_mile_pop', 'x7_min_pop', 'x1_2_mile_pop',
'dayofweek', 'quarter', 'month', 'year', 'dayofyear',
'dayofmonth', 'lag1', 'isFuture']
# Split the data into features (X) and target (y)
X = df[FEATURES]
y = df[TARGET]
# Create the XGBoost model
reg = xgb.XGBRegressor(base_score=0.5,
booster='gbtree',
n_estimators=500,
objective='reg:linear',
max_depth=3,
learning_rate=0.01)
# Fit the model
reg.fit(X, y, eval_set=[(X, y)], verbose=100)
Out[173]:
Creating Future Dataframe¶
In [174]:
future = pd.date_range('2023-08-17','2024-12-31', freq='D')
future_df = pd.DataFrame(index=future)
future_df['isFuture'] = True
df['isFuture'] = False
df_and_future = pd.concat([df, future_df])
df_and_future = create_features(df_and_future)
df_and_future = add_lags(df_and_future)
In [175]:
future_w_features = df_and_future.query('isFuture').copy()
In [176]:
future_w_features.head()
Out[176]:
In [177]:
future_w_features = future_w_features.drop('daily_yoy_ndt.total_inside_sales', axis=1)
In [178]:
future_w_features['pred'] = reg.predict(future_w_features)
In [179]:
future_w_features.index = pd.to_datetime(future_w_features.index) # Make sure index is datetime
plt.figure(figsize=(10, 5))
plt.plot(future_w_features.index, future_w_features['pred'], color='r', marker='o', linestyle='-', linewidth=1, markersize=3)
plt.title('Future Predictions For Inside Sales ')
plt.xlabel('Date')
plt.ylabel('Inside Sales')
plt.grid(True)
plt.show()
Combining Future Prediction with the Historical Data¶
In [180]:
plt.figure(figsize=(15, 5))
# Plot the entire historical data
encoded_df['daily_yoy_ndt.total_inside_sales'].plot(label='Historical Data', color='b')
plt.axvline('2023-08-01', color='black', ls='--')
# Plot future predictions
plt.plot(future_w_features.index, future_w_features['pred'], label='Future Predictions', color='red', linestyle='dashed')
# Set title and labels
plt.title('Historical Data and Future Predictions of Total Inside Sales')
plt.xlabel('Date')
plt.ylabel('Value')
# Add legend
plt.legend()
# Show plot
plt.show()
Displaying the data for future predictions¶
In [181]:
start_date = '2023-08-17'
end_date = pd.to_datetime(start_date) + pd.DateOffset(days=502)
date_range = pd.date_range(start=start_date, end=end_date, freq='D')
# Create a DataFrame for the forecasted sales
forecasted_sales_df = pd.DataFrame({'Date': date_range, 'Sales': future_w_features['pred']})
# Display the DataFrame
print(forecasted_sales_df)
C.Target Variable = Diesel¶
Plotting the graph for sale of diesel over the years¶
In [182]:
average_sales = encoded_df['diesel'].mean()
encoded_df['diesel'].plot(style=".", figsize=(15,5), title="Diesel Sale over the Years")
plt.axhline(y=average_sales, color='r', linestyle='-', label='Average Sales')
plt.legend()
plt.show()
Train Test Split¶
In [183]:
#splitting test as date less than 1st jan 2023 and train as later
train= encoded_df['diesel'].loc[encoded_df.index< '2023-01-01']
test= encoded_df['diesel'].loc[encoded_df.index>= '2023-01-01']
fig, ax= plt.subplots(figsize=(15,5))
train.plot(ax=ax, label='training set')
test.plot(ax=ax, label='Test set')
ax.axvline('2023-01-01',color='black',ls='--')
ax.legend(['Training Set','Test Set'])
plt.show()
Splitting data into training and testing sets¶
In [184]:
# Features (X) and Target (y)
X = encoded_df.drop(columns=['diesel'])
y = encoded_df['diesel']
# Split data into training and testing sets
train_end_date = '2023-01-01'
X_train = X[X.index < train_end_date]
X_test = X[X.index >= train_end_date]
y_train = y[y.index < train_end_date]
y_test = y[y.index >= train_end_date]
We now have X_train, X_test, y_train, and y_test ready for modeling¶
In [185]:
reg = xgb.XGBRegressor(base_score=0.5, booster='gbtree',
n_estimators=1000,
early_stopping_rounds=50,
objective='reg:linear',
max_depth=3,
learning_rate=0.01)
reg.fit(X_train, y_train,
eval_set=[(X_train, y_train), (X_test, y_test)],
verbose=100)
Out[185]:
Identifying Top 10 features impacting the Diesel Sale¶
In [186]:
fi = pd.DataFrame(data=reg.feature_importances_,
index=reg.feature_names_in_,
columns=['importance'])
top_10_features = fi.sort_values(by='importance', ascending=False).head(10)
top_10_features.plot(kind='barh', title='Top 10 Feature Importance')
Out[186]:
Plotting the graph to depict Train and Test Split¶
In [187]:
train = encoded_df[['diesel','hi_flow_lanes_fueling_positions','x1_mile_emp',
'womens_toilet_count','daily_yoy_ndt.total_food_service','rv_lanes_fueling_positions','site_id_msba','traditional_forecourt_fueling_positions','x1_mile_pop']].loc[encoded_df.index < '2023-01-01']
test = encoded_df[['diesel','hi_flow_lanes_fueling_positions','x1_mile_emp',
'womens_toilet_count','daily_yoy_ndt.total_food_service','rv_lanes_fueling_positions','site_id_msba','traditional_forecourt_fueling_positions','x1_mile_pop']].loc[encoded_df.index >= '2023-01-01']
fig, ax = plt.subplots(figsize=(15, 5))
train['diesel'].plot(ax=ax, label='Training Set', title='Data Train/Test Split')
test['diesel'].plot(ax=ax, label='Test Set')
ax.axvline('2023-01-01', color='black', ls='--')
ax.legend(['Training Set', 'Test Set'])
plt.show()
Time series cross validation¶
In [188]:
tss = TimeSeriesSplit(n_splits=5)
df = encoded_df[['diesel','hi_flow_lanes_fueling_positions','x1_mile_emp',
'womens_toilet_count','daily_yoy_ndt.total_food_service','rv_lanes_fueling_positions','site_id_msba','traditional_forecourt_fueling_positions','x1_mile_pop']].sort_index()
In [189]:
fig, axs = plt.subplots(5, 1, figsize=(15, 15), sharex=True)
fold = 0
for train_idx, val_idx in tss.split(df):
train = df.iloc[train_idx]
test = df.iloc[val_idx]
train['diesel'].plot(ax=axs[fold],
label='Training Set',
title=f'Data Train/Test Split Fold {fold}')
test['diesel'].plot(ax=axs[fold],
label='Test Set')
axs[fold].axvline(test.index.min(), color='black', ls='--')
fold += 1
plt.show()
Creating time series features based on time series index¶
In [190]:
def create_features(df):
df = df.copy()
df['dayofweek'] = df.index.dayofweek
df['quarter'] = df.index.quarter
df['month'] = df.index.month
df['year'] = df.index.year
df['dayofyear'] = df.index.dayofyear
df['dayofmonth'] = df.index.day
return df
df = create_features(df)
Defining Lags¶
In [191]:
def add_lags(df):
target_map= df['diesel'].to_dict()
df['lag1'] = (df.index - pd.Timedelta('364 days')).map(target_map)
return df
In [192]:
df = add_lags(df)
In [193]:
df.tail()
Out[193]:
Fitting the Model to calculate Evaluation Metrics¶
In [194]:
tss = TimeSeriesSplit(n_splits=5)
df = df.sort_index()
fold = 0
preds = []
scores = []
for train_idx, val_idx in tss.split(df):
train = df.iloc[train_idx]
test = df.iloc[val_idx]
TARGET = 'diesel'
X_train = train
y_train = train[TARGET]
X_test = test
y_test = test[TARGET]
reg = xgb.XGBRegressor(base_score=0.5, booster='gbtree',
n_estimators=1000,
early_stopping_rounds=50,
objective='reg:linear',
max_depth=3,
learning_rate=0.01)
reg.fit(X_train, y_train,
eval_set=[(X_train, y_train), (X_test, y_test)],
verbose=100)
y_pred = reg.predict(X_test)
preds.append(y_pred)
score = np.sqrt(mean_squared_error(y_test, y_pred))
scores.append(score)
Calculation of evaluation metrics¶
In [195]:
mae = mean_absolute_error(y_test, y_pred)
mse = mean_squared_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)
# Print the evaluation metrics
print(f'Mean Absolute Error (MAE): {mae:.2f}')
print(f'Mean Squared Error (MSE): {mse:.2f}')
print(f'R-squared (R2) Score: {r2:.2f}')
In [196]:
print(f'Score across folds {np.mean(scores):0.4f}')
print(f'Fold scores:{scores}')
Retraining all data¶
In [197]:
future = pd.date_range('2023-08-01','2024-12-31', freq='D')
future_df = pd.DataFrame(index=future)
df['isFuture'] = False
# Define your target variable
TARGET = 'diesel'
# Define features (excluding the target and 'isFuture' column)
FEATURES = ['hi_flow_lanes_fueling_positions','x1_mile_emp',
'womens_toilet_count','daily_yoy_ndt.total_food_service','rv_lanes_fueling_positions','site_id_msba','traditional_forecourt_fueling_positions','x1_mile_pop',
'dayofweek', 'quarter', 'month', 'year', 'dayofyear',
'dayofmonth','lag1', 'isFuture']
# Split the data into features (X) and target (y)
X = df[FEATURES]
y = df[TARGET]
# Create the XGBoost model
reg = xgb.XGBRegressor(base_score=0.5,
booster='gbtree',
n_estimators=500,
objective='reg:linear',
max_depth=3,
learning_rate=0.01)
# Fit the model
reg.fit(X, y, eval_set=[(X, y)], verbose=100)
Out[197]:
Creating Future Dataframe¶
In [198]:
future = pd.date_range('2023-08-17','2024-12-31', freq='D')
future_df = pd.DataFrame(index=future)
future_df['isFuture'] = True
df['isFuture'] = False
df_and_future = pd.concat([df, future_df])
df_and_future = create_features(df_and_future)
df_and_future = add_lags(df_and_future)
In [199]:
future_w_features = df_and_future.query('isFuture').copy()
In [200]:
future_w_features.head()
Out[200]:
In [201]:
future_w_features = future_w_features.drop('diesel', axis=1)
In [202]:
future_w_features['pred'] = reg.predict(future_w_features)
In [203]:
future_w_features.index = pd.to_datetime(future_w_features.index) # Make sure index is datetime
plt.figure(figsize=(10, 5))
plt.plot(future_w_features.index, future_w_features['pred'], color='r', marker='o', linestyle='-', linewidth=1, markersize=3)
plt.title('Future Predictions For Diesel ')
plt.xlabel('Date')
plt.ylabel('Diesel')
plt.grid(True)
plt.show()
Combining Future Prediction with the Historical Data¶
In [204]:
plt.figure(figsize=(15, 5))
encoded_df['diesel'].plot(label='Historical Data', color='b')
plt.axvline('2023-08-01', color='black', ls='--')
# Plot future predictions
plt.plot(future_w_features.index, future_w_features['pred'], label='Future Predictions', color='red', linestyle='dashed')
# Set title and labels
plt.title('Historical Data and Future Predictions of Diesel Sales')
plt.xlabel('Date')
plt.ylabel('Value')
plt.legend()
plt.show()
Displaying the data for future predictions¶
In [205]:
start_date = '2023-08-17'
end_date = pd.to_datetime(start_date) + pd.DateOffset(days=502)
date_range = pd.date_range(start=start_date, end=end_date, freq='D')
# Create a DataFrame for the forecasted sales
forecasted_sales_df = pd.DataFrame({'Date': date_range, 'Sales': future_w_features['pred']})
# Display the DataFrame
print(forecasted_sales_df)
D.Target Variable = Unleaded¶
Plotting the graph for sale of Unleaded over the years¶
In [206]:
average_sales = encoded_df['unleaded'].mean()
encoded_df['unleaded'].plot(style=".", figsize=(15,5), title="Unleaded Sale over the Years")
plt.axhline(y=average_sales, color='r', linestyle='-', label='Average Sales')
plt.legend()
plt.show()
Train Test Split¶
In [207]:
train= encoded_df['unleaded'].loc[encoded_df.index< '2023-01-01']
test= encoded_df['unleaded'].loc[encoded_df.index>= '2023-01-01']
fig, ax= plt.subplots(figsize=(15,5))
train.plot(ax=ax, label='training set')
test.plot(ax=ax, label='Test set')
ax.axvline('2023-01-01',color='black',ls='--')
ax.legend(['Training Set','Test Set'])
plt.show()
Splitting data into training and testing sets¶
In [208]:
X = encoded_df.drop(columns=['unleaded'])
y = encoded_df['unleaded']
# Split data into training and testing sets
train_end_date = '2023-01-01'
X_train = X[X.index < train_end_date]
X_test = X[X.index >= train_end_date]
y_train = y[y.index < train_end_date]
y_test = y[y.index >= train_end_date]
We now have X_train, X_test, y_train, and y_test ready for modeling¶
In [209]:
reg = xgb.XGBRegressor(base_score=0.5, booster='gbtree',
n_estimators=1000,
early_stopping_rounds=50,
objective='reg:linear',
max_depth=3,
learning_rate=0.01)
reg.fit(X_train, y_train,
eval_set=[(X_train, y_train), (X_test, y_test)],
verbose=100)
Out[209]:
Identifying Top 10 features impacting the Food Sale¶
In [210]:
fi = pd.DataFrame(data=reg.feature_importances_,
index=reg.feature_names_in_,
columns=['importance'])
top_10_features = fi.sort_values(by='importance', ascending=False).head(10)
top_10_features.plot(kind='barh', title='Top 10 Feature Importance')
Out[210]:
Plotting the graph to depict Train and Test Split¶
In [211]:
train = encoded_df[['unleaded','mens_toilet_count','x1_mile_income', 'traditional_forecourt_stack_type','square_feet','lottery','years_since_last_project','x7_min_emp','total_sales']].loc[encoded_df.index < '2023-01-01']
test = encoded_df[['unleaded','mens_toilet_count','x1_mile_income', 'traditional_forecourt_stack_type','square_feet','lottery','years_since_last_project','x7_min_emp','total_sales']].loc[encoded_df.index >= '2023-01-01']
fig, ax = plt.subplots(figsize=(15, 5))
train['unleaded'].plot(ax=ax, label='Training Set', title='Data Train/Test Split')
test['unleaded'].plot(ax=ax, label='Test Set')
ax.axvline('2023-01-01', color='black', ls='--')
ax.legend(['Training Set', 'Test Set'])
plt.show()
Time series cross validation¶
In [212]:
tss = TimeSeriesSplit(n_splits=5)
df = encoded_df[['unleaded','mens_toilet_count','x1_mile_income', 'traditional_forecourt_stack_type','square_feet','lottery','years_since_last_project','x7_min_emp','total_sales']].sort_index()
In [213]:
fig, axs = plt.subplots(5, 1, figsize=(15, 15), sharex=True)
fold = 0
for train_idx, val_idx in tss.split(df):
train = df.iloc[train_idx]
test = df.iloc[val_idx]
train['unleaded'].plot(ax=axs[fold],
label='Training Set',
title=f'Data Train/Test Split Fold {fold}')
test['unleaded'].plot(ax=axs[fold],
label='Test Set')
axs[fold].axvline(test.index.min(), color='black', ls='--')
fold += 1
plt.show()
Creating time series features based on time series index¶
In [214]:
def create_features(df):
df = df.copy()
df['dayofweek'] = df.index.dayofweek
df['quarter'] = df.index.quarter
df['month'] = df.index.month
df['year'] = df.index.year
df['dayofyear'] = df.index.dayofyear
df['dayofmonth'] = df.index.day
return df
df = create_features(df)
Defining Lags¶
In [215]:
def add_lags(df):
target_map= df['unleaded'].to_dict()
df['lag1'] = (df.index - pd.Timedelta('364 days')).map(target_map)
return df
In [216]:
df = add_lags(df)
In [217]:
df.tail()
Out[217]:
Fitting the Model to calculate Evaluation Metrics¶
In [218]:
tss = TimeSeriesSplit(n_splits=5)
df = df.sort_index()
fold = 0
preds = []
scores = []
for train_idx, val_idx in tss.split(df):
train = df.iloc[train_idx]
test = df.iloc[val_idx]
TARGET = 'unleaded'
X_train = train
y_train = train[TARGET]
X_test = test
y_test = test[TARGET]
reg = xgb.XGBRegressor(base_score=0.5, booster='gbtree',
n_estimators=1000,
early_stopping_rounds=50,
objective='reg:linear',
max_depth=3,
learning_rate=0.01)
reg.fit(X_train, y_train,
eval_set=[(X_train, y_train), (X_test, y_test)],
verbose=100)
y_pred = reg.predict(X_test)
preds.append(y_pred)
score = np.sqrt(mean_squared_error(y_test, y_pred))
scores.append(score)
Calculation of evaluation metrics¶
In [219]:
mae = mean_absolute_error(y_test, y_pred)
mse = mean_squared_error(y_test, y_pred)
r2 = r2_score(y_test, y_pred)
print(f'Mean Absolute Error (MAE): {mae:.2f}')
print(f'Mean Squared Error (MSE): {mse:.2f}')
print(f'R-squared (R2) Score: {r2:.2f}')
In [220]:
print(f'Score across folds {np.mean(scores):0.4f}')
print(f'Fold scores:{scores}')
Retraining all data¶
In [221]:
future = pd.date_range('2023-08-01','2024-12-31', freq='D')
future_df = pd.DataFrame(index=future)
df['isFuture'] = False
TARGET = 'unleaded'
FEATURES = ['mens_toilet_count','x1_mile_income', 'traditional_forecourt_stack_type','square_feet','lottery','years_since_last_project','x7_min_emp','total_sales',
'dayofweek', 'quarter', 'month', 'year', 'dayofyear',
'dayofmonth', 'lag1', 'isFuture']
X = df[FEATURES]
y = df[TARGET]
reg = xgb.XGBRegressor(base_score=0.5,
booster='gbtree',
n_estimators=500,
objective='reg:linear',
max_depth=3,
learning_rate=0.01)
reg.fit(X, y, eval_set=[(X, y)], verbose=100)
Out[221]:
Creating Future Dataframe¶
In [222]:
future = pd.date_range('2023-08-17','2024-12-31', freq='D')
future_df = pd.DataFrame(index=future)
future_df['isFuture'] = True
df['isFuture'] = False
df_and_future = pd.concat([df, future_df])
df_and_future = create_features(df_and_future)
df_and_future = add_lags(df_and_future)
In [223]:
future_w_features = df_and_future.query('isFuture').copy()
In [224]:
future_w_features.head()
Out[224]:
In [225]:
future_w_features = future_w_features.drop('unleaded', axis=1)
In [226]:
future_w_features['pred'] = reg.predict(future_w_features)
In [227]:
future_w_features.index = pd.to_datetime(future_w_features.index) # Make sure index is datetime
plt.figure(figsize=(10, 5))
plt.plot(future_w_features.index, future_w_features['pred'], color='r', marker='o', linestyle='-', linewidth=1, markersize=3)
plt.title('Future Predictions For Unleaded ')
plt.xlabel('Date')
plt.ylabel('Unleaded')
plt.grid(True)
plt.show()
Combining Future Prediction with the Historical Data¶
In [228]:
plt.figure(figsize=(15, 5))
encoded_df['unleaded'].plot(label='Historical Data', color='b')
plt.axvline('2023-08-17', color='black', ls='--')
plt.plot(future_w_features.index, future_w_features['pred'], label='Future Predictions', color='red', linestyle='dashed')
plt.title('Historical Data and Future Predictions of Unleaded Sales')
plt.xlabel('Date')
plt.ylabel('Value')
plt.legend()
plt.show()
Displaying the data for future predictions¶
In [229]:
start_date = '2023-08-17'
end_date = pd.to_datetime(start_date) + pd.DateOffset(days=502)
date_range = pd.date_range(start=start_date, end=end_date, freq='D')
forecasted_sales_df = pd.DataFrame({'Date': date_range, 'Sales': future_w_features['pred']})
print(forecasted_sales_df)
II. PROPHET¶
In [230]:
#!pip install pystan==2.19.1.1
In [231]:
#pip install --upgrade numpy
In [232]:
from prophet import Prophet
In [233]:
df_merged['calendar.calendar_day_date'] = t_data.index
df_merged.set_index('calendar.calendar_day_date', inplace=True)
In [235]:
import pandas as pd
categorical_columns = df_merged.select_dtypes(include=['object']).columns
encoded_df = pd.get_dummies(df_merged, columns=categorical_columns)
In [236]:
encoded_df.head()
Out[236]:
In [237]:
split_date = '2023-01-01'
inside_train = encoded_df['daily_yoy_ndt.total_inside_sales'].loc[encoded_df.index <= split_date].copy()
inside_test = encoded_df['daily_yoy_ndt.total_inside_sales'].loc[encoded_df.index > split_date].copy()
# Plot train and test so you can see where we have split
(inside_test
.rename_axis('Date')
.rename('TEST SET')
.to_frame()
.join(inside_train.rename_axis('Date').rename('TRAINING SET').to_frame(),
how='outer')
.plot(figsize=(10, 5), title='Inside Sales', style='.', ms=2)
)
plt.show()
In [238]:
inside_train_prophet= inside_train.reset_index().rename(columns={'calendar.calendar_day_date':'ds','daily_yoy_ndt.total_inside_sales':'y'})
In [239]:
inside_train_prophet.head()
Out[239]:
In [240]:
model= Prophet()
model.fit(inside_train_prophet)
Out[240]:
In [241]:
inside_test_prophet= inside_test.reset_index().rename(columns={'calendar.calendar_day_date':'ds','daily_yoy_ndt.total_inside_sales':'y'})
inside_sales_tst_forecast= model.predict(inside_test_prophet)
In [242]:
inside_sales_tst_forecast.head()
Out[242]:
In [243]:
fig, ax = plt.subplots(figsize=(10, 5))
fig = model.plot(inside_sales_tst_forecast, ax=ax)
plt.show()
In [244]:
fig = model.plot_components(inside_sales_tst_forecast)
plt.show()
In [245]:
f, ax = plt.subplots(figsize=(15,5))
ax.scatter(inside_test.index, inside_test, color='#FFB6C1')
fig = model.plot(inside_sales_tst_forecast, ax=ax)
In [246]:
f, ax = plt.subplots(figsize=(15,5))
ax.scatter(inside_test.index, inside_test, color='r')
fig = model.plot(inside_sales_tst_forecast, ax=ax)
ax.set_xlim(left=pd.to_datetime('2023-01-01'), right=pd.to_datetime('2023-01-31'))
ax.set_ylim(bottom=0, top=6000)
plt.suptitle('January 2023 Forecast vs Actual')
plt.show()
In [247]:
#Jan 2023 first week forecasts in detail
f, ax = plt.subplots(figsize=(15,5))
ax.scatter(inside_test.index, inside_test, color='r')
fig = model.plot(inside_sales_tst_forecast, ax=ax)
ax.set_xlim(left=pd.to_datetime('2023-01-01'), right=pd.to_datetime('2023-01-07'))
ax.set_ylim(bottom=0, top=6000)
plt.suptitle('January 2023 Forecast vs Actual')
plt.show()
In [248]:
#Evaluate the model with error metrics
from sklearn.metrics import mean_squared_error
np.sqrt(mean_squared_error(y_true=inside_test,y_pred=inside_sales_tst_forecast['yhat']))
Out[248]:
In [249]:
def mean_absolute_percentage_error(y_true, y_pred):
y_true, y_pred = np.array(y_true), np.array(y_pred)
return np.mean(np.abs((y_true - y_pred) / y_true)) * 100
mean_absolute_percentage_error(y_true=inside_test, y_pred=inside_sales_tst_forecast['yhat'])
Out[249]:
In [250]:
future= model.make_future_dataframe(periods=730)
forecast=model.predict(future)
In [251]:
fig = model.plot(forecast, figsize=(15, 5))
plt.title('Prophet Forecast')
plt.xlabel('Date')
plt.ylabel('Value')
plt.show()
In [252]:
forecast[['ds','yhat']].tail(100)
Out[252]:
In our analysis using the Prophet model, we've noticed that the one-year forecast for the Inside Sales target variable exhibits a decreasing trend, which contradicts the actual values represented in light red, indicating a potential overfitting issue. This discrepancy between predicted and actual values can be attributed to certain limitations of the Prophet model that became evident during our dataset testing.¶
Some of the notable limitations of the Prophet model include its simplicity in handling complex interactions and dependencies within the time series data. The model's approach may not effectively capture intricate relationships among various influencing factors, leading to suboptimal predictions.¶
Additionally, Prophet makes certain assumptions about seasonality, primarily focusing on weekly and yearly patterns. However, if the actual seasonal patterns in the data follow a different trajectory, the model may struggle to accurately account for these variations.¶
4. Model Performance¶
XGBoost Model Performance Food Sales Mean Absolute Error (MAE): 2.51 Mean Squared Error (MSE): 34.41
Inside Sales Mean Absolute Error (MAE): 6.68 Mean Squared Error (MSE): 423.45
Diesel: Mean Absolute Error (MAE): 12.91 Mean Squared Error (MSE): 2676.26
Unleaded: Mean Absolute Error (MAE): 5.83 Mean Squared Error (MSE): 166.44
The XGBoost model's performance across the four different target variables, namely Food Sales, Inside Sales, Diesel, and Unleaded, reveals interesting insights. Firstly, when it comes to the Mean Absolute Error (MAE), we observe that the Food Sales and Unleaded categories have the lowest errors, with MAE values of 2.51 and 5.83, respectively. This indicates that the model's predictions for these two variables are relatively close to the actual values. In contrast, the Inside Sales and Diesel categories have higher MAE values of 6.68 and 12.91, respectively, suggesting more significant prediction errors.¶
Additionally, the Mean Squared Error (MSE) provides a measure of the model's predictive accuracy, where lower values indicate better performance. In this context, Food Sales and Unleaded again stand out with the lowest MSE values of 34.41 and 166.44, respectively. Conversely, Inside Sales and Diesel have considerably higher MSE values of 423.45 and 2676.26, implying that the model's predictions for these variables deviate more from the actual values.¶
Given these results, XGBoost appears to be best suited for predicting Food Sales and Unleaded fuel sales. These two variables exhibit the lowest MAE and MSE, indicating the model's better performance and accuracy in predicting them. However, when it comes to Inside Sales and Diesel, the higher MAE and MSE values suggest that the model may not perform as well for these variables, and alternative modeling approaches may be more appropriate. It's essential to consider the specific context and business goals to determine which target variable is of higher importance and, consequently, which one the XGBoost model is best suited for. The XGBoost model delivered strong performance during training, but when applied to the test set, we observed some discrepancies in the predictions for Inside Sales and Diesel. Consequently, we are actively exploring alternative models that can offer improved predictive accuracy for these specific target variables.¶
5. Results¶
In our analysis, we've employed both the Prophet and XGBoost models to forecast various target variables, shedding light on their respective strengths and limitations. The Prophet model exhibited a notable drawback in its one-year forecast for the Inside Sales target variable, where it portrayed a decreasing trend that contradicted the actual values, indicating a potential overfitting issue. This discrepancy can be attributed to the model's inherent limitations, particularly its simplistic approach in capturing complex interactions and dependencies within time series data. Additionally, the model's assumptions about seasonality, primarily focusing on weekly and yearly patterns, may lead to inaccuracies when the actual data follows a different seasonal trajectory.¶
Conversely, the XGBoost model's performance across four different target variables, namely Food Sales, Inside Sales, Diesel, and Unleaded, showcased intriguing results. It revealed that the Food Sales and Unleaded categories had the lowest Mean Absolute Error (MAE) and Mean Squared Error (MSE) values, indicating the model's ability to predict these variables with high accuracy. In contrast, Inside Sales and Diesel displayed higher MAE and MSE values, implying a more significant divergence between predicted and actual values.¶
While XGBoost demonstrates improved performance compared to Prophet, it's essential to note that the high R2 value suggests potential overfitting in the model. We are currently exploring ways to fine-tune the model to achieve more optimal and reliable results.¶
6. Contributions¶
Krishna Chaitanya: XgBoost Model(Model built for Inside Sales), Correlation Matrix, Seasonality, Analysis of unleaded and diesel sales, Results
Anukriti Raj: Arimax Model, XgBoost Model(Model built for Food Sales) Feature Engineering ,Seasonality, Correlation Matrix,Analysis of unleaded and diesel sales, Results
Shivi Shrivastav: Introduction, XgBoost Model(Model built for Unleaded Sales) Table of Contents, Seasonality, Analysis of in-store and food sales sales, Results
Litzy Carbajal: Introduction, XgBoost Model(Model built for Diesel Sales) ,Seasonality, Analysis of in-store and food sales, Results